%% This code estimates risk premia for each commodity using a K-factor model.
%% For each factor, the model includes one jump component without any jump size restriction, i.e., over the full range [-inf, inf].
%% The estimates are saved in 'CM_J_Model_', num2str(K), '.mat', which are later used to generate tables and figures.
%%
%% The same code can also estimate a model that does not separate jump risk premia from continuous risk,
%% by setting the JumpIndicator flag to false.
%%
%%
%% 'CM_factors.mat' contains the constructed factor returns at hourly frequency.
%% - Column 1: Date
%% - Column 2: Time index, where 0 indicates the return from 11pm to midnight local time.
%%             Prices are forward-filled using the last available price when not observed.
%% - Columns 3–7: Returns for the five constructed factors.
%% - Column 8: Synthetic day index, constructed with 24 observations per day.
%% - Column 9: Synthetic hour index (0 to 23).
%%
%% 'CM_returns_sample.mat' contains a small sample of 22 commodity futures returns used in the empirical analysis (Nov 2015).
%% The full dataset is restricted and cannot be shared. The data are orgnized in the same format as 'CM_factors.mat'.
%%
%% All estimates are based on synthetic days to mitigate issues arising from irregular trading hours, holidays, or changes in session length.

clear;

addpath('../FX/Functions') % CM and FX share the same functions

load('CM_factors.mat');
load('CM_returns.mat');
%% Continuous Factor and Jump Specification


Models = {{'3'},{'3','2'},{'3','2','5'},{'3','2','5','1'},{'3','2','5','1','4'}};

JumpIndicator = true;%false;


for model    = 1:length(Models)  % 
    disp(factors_name(str2double(Models{model})));

K        =  length(Models{model});

if JumpIndicator

    B_kl = repmat([-inf,inf],K,1);
    L_k  =  ones(K,1);

else

    B_kl=[];
    L_k=zeros(K,1);
    
end
    H    =   sum(L_k);

%% Empirical Setting
numofdays = rets(end,end-1);
nobsperday = (max(factors(:,end))+1);
numoftickers = size(rets,2) - 4;

%% 

qn       =   21*nobsperday;
Delta_n  =   1/252/nobsperday;
M        =   numoftickers;
numofwindows = floor(numofdays/21);
n = numofwindows * 21 * nobsperday;
rets = rets(1:n,:);
factors =factors(1:n,:);
Tn       =  n*Delta_n;
 

lambdat_hat = zeros(numofwindows, K+H);
Vnt_hat     = zeros(K+H, K+H, numofwindows);
betat_hat   = zeros(M, K+H, numofwindows);
yt_hat      = zeros(M, numofwindows);
R2_t        = zeros(numofwindows,1);

dPt    =  rets(:,3:(end-2))';
dFt    =  factors(:,2+str2double(Models{model}))';

[beta_J_hat, R2_j] = Jump_Beta(dFt, dPt, n, Delta_n, nobsperday, B_kl, L_k);

for i = 1:numofwindows-1
    disp(i)
    dPt_j    =  rets((qn*(i-1)+1):(qn*i),3:(end-2))'; % Current-Period return
    dFt_j    =  factors((qn*(i-1)+1):(qn*i),2+str2double(Models{model}))';
    dPt_jp1    =  rets((qn*(i)+1):(qn*(i+1)),3:(end-2))'; % Next-Period return 
   [lambdat_hat(i,:), Vnt_hat(:,:,i), betat_hat(:,:,i), yt_hat(:,i), R2_t(i,:)] = AJX_HF(dFt_j, dPt_j, dPt_jp1, beta_J_hat, qn, Delta_n, nobsperday, Tn);
end

Lambda_hat_AJX = sum(lambdat_hat,1)/(floor(n/qn)*qn*Delta_n);
Vn_AJX         = sum(Vnt_hat,3)/(floor(n/qn)*qn*Delta_n)^2;

if JumpIndicator
    save(['CM_J_Model_',num2str(model),'.mat'],'lambdat_hat', 'Lambda_hat_AJX', 'Vnt_hat', 'Vn_AJX', 'betat_hat', 'yt_hat', 'R2_t', 'qn', 'K', 'H', 'Delta_n')
else
    save(['CM_Model_',num2str(model),'.mat'],'lambdat_hat', 'Lambda_hat_AJX', 'Vnt_hat', 'Vn_AJX', 'betat_hat', 'yt_hat', 'R2_t', 'qn', 'K', 'H', 'Delta_n')
end

end